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Abstract 

Assimilation of surface geomagnetic observations and geodynamo models has 
advanced very quickly in recent years. However, compared to advanced data 
assimilation systems in meteorology, geomagnetic data assimilation (GDAS) is 
still in an early stage. Among many challenges ranging from data to models is the 
disparity between the short observation records and the long time scales of the 
core dynamics. To better utilize available observational information, we have 
made an effort in this study to directly assimilate the Gauss coefficients of both 
the core field and its secular variation (SV) obtained via global geomagnetic field 
modeling, aiming at understanding the dynamical responses of the core fluid to 
these additional observational constraints. Our studies show that the SV 
assimilation helps significantly to shorten the dynamo model spin-up process. The 
flow beneath the core-mantle boundary (CMB) responds significantly to the 
observed field and its SV. The strongest responses occur in the relatively small 
scale flow (of the degrees L « 30 in spherical harmonic expansions). This part of 
the flow includes the axisymmetric toroidal flow (of order m = 0) and 
non-axisymmetric poloidal flow with m > 5. These responses can be used to 
better understand the core flow and, in particular, to improve accuracies of 
predicting geomagnetic variability in future. 
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Introduction 

Geomagnetic field observed at the Earth’s surface varies significantly in time: its 
temporal scales range from minutes to geological time scales. Though it was first 
noticed by mankind over 5000 years ago (Roberts, 1992), and its origin was sought 
as early as 800 years ago (Dibner Library 1980), the modern theory that the geo- 
magnetic held is generated and maintained by convective how in the Earth’s outer 
core (geodynamo) was originated from the seminal work of Larmor (1919). Success- 
ful numerical simulation of the geodynamo was hrst carried out by Glatzmaier and 
Roberts (1995), and then followed by Kageyama and Sato (1997), and by Kuang 
and Bloxham (1997). Christensen et al (2010) provided a comprehensive summary 
of numerical geodynamo solutions and their relevances to geomagnetic observations. 

Assimilation of geomagnetic observations with numerical geodynamo models 
started less than a decade ago. Sun et al (2007) and Fournier et al (2007) used 
simplified magnetohydrodynamic (MHD) systems and synthetic data tested the 
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applicability of assimilation of sparse magnetic data. Liu et al (2007) first used ob- 
servation system simulation experiments (OSSEs) with a full dynamo and demon- 
strated clearly that one could use assimilation of magnetic field at the surface to 
estimate the dynamo state deep in the fluid core. Kuang et al (2008) published 
the first working geomagnetic data assimilation system MoSST_DAS in which the 
Gauss coefficients of various geomagnetic and paleomagnetic field models are as- 
similated with their MoSST geodynamo model (Kuang and Chao, 2003; Jiang and 
Kuang, 2008) for estimation of the core state and prediction of geomagnetic field. 
Kuang et al (2009) then used this assimilation system and 100 years of the Gauss 
coefficients from GUFM1 (Jackson et al 2000) and CM4 (Sabaka et al 2004) to un- 
derstand the responses of the core state to surface geomagnetic observations, and 
their implications to core state estimation and SV prediction. We refer the reader to 
Fournier et al (2010) for a comprehensive review of the data assimilation algorithms 
for geomagnetic data assimilation (GDAS) and some of the early results. 

Rapid advances have occurred in multiple facets of GDAS. Several independent 
assimilation systems have been developed to understand better the core dynamical 
state. For example, Aubert and Fournier (2011), and Fournier et al (2011, 2013) 
carried out OSSEs with synthetic observations and numerical dynamo models to 
examine possibilities of core state determination. Aubert (2013, 2014) investigated 
possibilities of inverting core state properties using the observed field and SV. In 
addition to the sequential data assimilation systems mentioned above, there are 
also efforts in developing GDAS systems based on variational data assimilation 
techniques. For example, Li et al (2011, 2014) have been continuing their effort on a 
new combined system of forward and adjoint systems. Encompassed application is 
the contributions of assimilation results to international geomagnetic reference field 
(IGRF) (Kuang et al, 2010), and efforts to determine field model error statistics 
(Gillet et al, 2013). 

Despite these advances, GDAS is still in an early stage similar to that of early 
numerical weather prediction (NWP) (for a more comprehensive review, see, e.g. 
Kalnay 2003). Many important questions are still to be fully answered, such as 
comprehensive assessment of numerical dynamo system biases, observation and core 
state covariances and error statistics, and the dynamic responses of dynamo state 
to the observed geomagnetic field. The latter is of in particular importance to the 
spin-up processes of the numerical models which, in turn, determine how fast and 
how close the numerical solutions can be pulled to the true state of the core. 

Concerns on the spin-up of the numerical models can be examined from the time 
scales of the observed field and of the numerical models. Global field model results 
from the past 400 years of geomagnetic data (e.g. Jackson et al, 2000; Sabaka et 
al, 2004, 2015; Olsen et al, 2006, 2014) show that the typical time scales T; of the 
degree l components (Stacy 1992; Hulot and Le Mouel 1994; Olsen et al 2006) 
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varies from over 1000 years for the dipole (l = 1) to less than 100 years for higher 
degrees (see Figure 1). In (1), (g ; m , h™) are the Gauss coefficients of the field, and 
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(< 7 ; m , Zij 71 ) are their first order time derivatives, i.e. the Gauss coefficients of the SV. 
Currently, the longest record for low degree (Z < 5) field coefficients is from the 
paleo/archeo magnetic data (e.g. Korte et al , 2011; Nilsson et al , 2014). The high 
quality coefficients for up to degree l < 8 could be obtained from historical and 
observatory data (Jackson et al, 2000). Very high quality coefficients for degrees 
l < 13 are obtained in the past 50 years with satellite magnetic data (Sabaka et 
al, 2004, 2015; Olsen et al, 2006, 2014). In summary, the data record is no more 
than 10 times of the typical time scales of the geomagnetic field. This brings the 
very concern on whether the observational record is sufficient to spin up numerical 
dynamo models. The model spin-up also has direct consequence on estimation of 
the core state. 

How could we improve geomagnetic data assimilation systems within the observa- 
tional limit? There are several areas for improvements. For example, improvements 
in global geomagnetic field modeling are needed since the Gauss coefficients from 
various field models have been used in most of the previous GDAS studies. Cur- 
rently there are many field models covering different epochs (e.g. Jackson et al, 
2000; Korte et al, 2011; Gillet et al 2013; Olsen et al, 2014; Sabaka et al, 2015). A 
unified field model covering the longest possible period could certainly reconcile dif- 
ferences in these models, and thus help greatly GDAS systems. There is an ongoing 
effort on constructing a unified global field model of the past millennium (private 
communication with Korte). The field model error statistical information of such 
unified field models, such as those in Gillet et al (2013), is also necessary for GDAS. 

Improvement in the assimilation algorithms could also help data utilization. Some 
efforts were made by Kuang et al (2010) in which a subset of the Gauss coefficients 
(of lower degrees) with much longer records are assimilated first to speed up the 
model, followed by assimilating those of higher degrees for the past 100 years. 
Tangborn and Kuang (2015) showed, via a set of experiments, that such assimilation 
methodology can have positive impact on core state, and improve accuracies of 
predicting the subset of the Gauss coefficients not assimilated. Another example 
is employment of ensemble Kalman Filtering (EnKF) approach (Evensen, 1994). 
Fournier et al (2011, 2013) used OSSEs to show the potential to speed up the 
transfer of information from geomagnetic data to the core state. But such speedy 
transfer depends on model errors (that are in general very large due to limitations of 
numerical dynamo models) not considered in their studies. It should also point out 
that GDAS is computationally very expensive. Such expense needs to be considered 
in the algorithm improvement. 

Another improvement is on exploiting and utilizing further geodynamic informa- 
tion embedded in surface geomagnetic measurements. An immediate candidate for 
such exploitation is the geomagnetic secular variation (SV), described by the first 
order time derivative (g™, h ™) of the Gauss coefficients since, as we will describe in 
the next section, they provide additional constraints on the core flow beneath the 
CMB, and on the radial variation of the magnetic field. The former is not new, as 
there is a long history of, started from Roberts and Scott (1965), core flow inver- 
sion from observed SV at the Earth’s surface via the “frozen-flux” approximation 
(in which the Ohmic dissipation beneath the CMB is ignored). However, this ap- 
proximation comes with the price: the core flow cannot be uniquely inverted (e.g. 
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Roberts and Scott 1965; Backus 1968). Thus additional constraints on core flow 
properties are necessary in such core flow inversion studies (for more complete re- 
views, please read, e.g., Holme, 2007; Kuang and Tangborn, 2011). If the Ohmic 
dissipation is retained (no “frozen flux” approximation), then the observed SV im- 
poses the constraints on the radial variation of the field in the core, as the latter is 
part of the magnetic induction. Since both field advection and Ohmic dissipation 
are included in geodynamo modeling, both kinds of constraints can be examined in 
MoSST_DAS or any other GDAS system without mathematical difficulties. 

Therefore, a natural expansion of data utilization in GDAS is to assimilate both 
the field and its SV, so that the embedded geodynamic constraints can be used to 
make more optimal analysis, thus speeding up the transport of information from 
the surface geomagnetic observations to the dynamical state in the outer core. Since 
the SV is not included in the state vector of numerical geodynamo models, it will 
be connected through a non-linear observation operator, T-L, which transforms the 
model state space to the observations space. Obviously H will depend on, among 
others, fundamental physical properties of the magnetic field. 

It should be pointed out here that assimilating the rate of change of geodynamic 
observables has been routinely used in numerical weather prediction (NWP). For 
example, precipitation rate, measured from a variety of satellite instruments is as- 
similated, despite not being a state variable in a GCM (Hou et al , 2000). It should 
also be pointed out that, in addition to core flow inversion (Roberts and Scott, 
1965), there are also attempts to invert core dynamical state with both the surface 
observations and the dynamo models (Aubert, 2013, 2014). The latter will benefit 
the SV assimilation. 

In this paper, we describe in detail the results from our recent effort on assimila- 
tion of both the field and its SV. These results, from a series of experiments, will 
demonstrate the improvement in prediction, and knowledge on core flow responses 
to the SV assimilation. The results also provide valuable information for further 
development in this direction. 

This paper is organized as follows: the numerical model details and the mathe- 
matical formulation for SV assimilation will be given in the next section. Followed 
are the experimental results we have with this assimilation approach. Discussions 
and plans for further improvements are presented in the last Section. 

Mathematical Description 

The mathematical formulation for SV assimilation depends on the numerical geody- 
namo models and the assimilation algorithms, in addition to the physics controlling 
the time variation of the magnetic field. In this section, we provide the mathemati- 
cal methodologies used in MoSST_DAS employed in this study (Kuang et al, 2008; 
Sun and Kuang, 2015). But, with some modifications, they can be applied to other 
GDAS systems. 

Dynamo State Vector and Geomagnetic Observation 

MoSST_DAS utilizes the MoSST core dynamics model for time integration of the 
magnetic field (Kuang et al, 2008; Sun and Kuang, 2015). In this system, the state 
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vector x 

X = (v,B,<5p) T 


(2) 


includes the velocity field v and the density anomaly Sg in the outer core r* < r < r c 
(ri and r c are the mean radii of the ICB and CMB, respectively); and the magnetic 
field B in the outer core, the electrically conducting inner core r < r* and the D”- 
layer r c < r < rd {^d is the mean radius at the top of the layer). The superscript 
“T” in (2) implies the transpose. The solid mantle above the D”-layer < r < r s 
(r s is the mean radius of the Earth’s surface) is electrically insulating. The whole 
system is defined in the reference frame fixed with the solid mantle. 

The velocity field v and the magnetic field B are decomposed into the poloidal 
and toroidal components, with the scalars described via spherical expansions 


(v, B) T = V x (T v ,T b f 1 T 


+ V x V x 


( Pv,Pb ) lr 


(p v ,T v ,p b ,T b ,6ef = J2 (vr^r^rjr^rfyr^^+c.c., 


0 <ra<Z 


( 3 ) 

( 4 ) 


where l r is the unit radial vector, 0 is the co-latitude, <f> is the longitude, are 
the fully normalized spherical harmonic functions of degree l and order m, Lm is 
the truncation order, and C.C. implies the complex conjugate part. P and T in (3) 
are called the poloidal and toroidal scalars. It is therefore convenient to write 


x 


= (x„ 


( 5 ) 


where the subsets are defined with the relevant spectral coefficients in (4), e.g., 


x b = {6™(r fc ) | 0 < r k < r d ; 0 < m < l < L M } T 


( 6 ) 


for the poloidal magnetic field. (5) and (6) can be different for other dynamo models. 

In geomagnetic field modeling, geomagnetic measurements are used to obtain the 
magnetic field B° originated from the core (simply called the geomagnetic field 
hereafter) that is described as 


B° = -Vtf , (7) 

4/ = r s ^ (— J (g™ cosrruj) + h™ sin m<j>) P™ (9) (8) 

0<m</ 

where P ; m is the Schmidt normalized associate Legendre polynomial of degree l and 
order to, {g r { 1 ■ h] n ) are the Gauss coefficients (slightly different from the standard 
notation), and L a is the maximum degree ( L a < 13 in general). Since these Gauss 
coefficients (g ; m , h are provided by different field models over the past 10000 
years (e.g. Jackson et al , 2000; Korte et al , 2005, 2011; Gillet et al , 2013; Olsen et 
al, 2014; Sabaka et al, 2015), they are used as the “observations” in our study. 
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By (3), (4), (7) and (8), we can obtain the relationship between (g™, h ™) in (8) 
and V" 1 in (4) via the radial component B r of the magnetic field B 


b? = -— = 


f)\ Xj J / r \ £+2 

E V + 1 Kt) (^ m cosm0- 


■hY 1 sin m<j>) PJ n (6) 


L 


0<m<Z 

La 


= "^= E 


l l± 1 l b Y<°) Yl m + c.c. 


( 9 ) 


0 <m<Z 


With the definitions of Y \ m and P™, (9) requires that 


27t(1 + SmO ) 


21 + 1 


1/2 


( 10 ) 


b? ( °\r) = r f(^)G m (gr-ihr), G m 
for r d < r < r s . The spectral coefficients of the SV are the time derivatives of (10): 
^ 0) (r)= r j(j) G m (gr- l hr) for r d <r<r s , (11) 


where (') means the time derivative. 


SV and Core State 

Geomagnetic observations only provide the time series of {gf 1 , h,™). The SV co- 
efficients (g™, h™) are actually derived. Assimilation of the SV thus raises two 
major concerns: could the SV be approximated as “instantaneously” measured, 
and whether it is redundant to the assimilation of the field? 

Answers to the first concern depend on the significance of numerical errors in SV 
calculation. Consider, for example, a central difference scheme is used, 

■ m(A _ 9f(* + St ) ~ 9?{t - 
91 W ~ 2 It ■ 

Then the relative numerical error is of order 


= O 


( To/tiY 


where r Q is the typical time intervals of data series, and p, defined in (1), is the 
typical time scales of the observed geomagnetic field. In general, t„ < 1 month in the 
field models using modern observatory and satellite data (e.g. Sabaka et al , 2004, 
2015; Olsen et al, 2006, 2014), while t; < 70 years (see Figure 1). Thus e n w 10~ 6 , 
which leads to an order 10~ 4 nT/year error in SV. On the other hand, the external 
field is several tens of nT at the Earth’s surface (Sabaka et al, 2015), and changes on 
the solar cycle (~ 11 years) and shorter time scales. Thus, e n is negligible compared 
to those arising from, e.g. separation of the external and the internal magnetic 
signals. One could then argue that both the field and its SV are “concurrently” 
measured. 

The redundancy is not an issue because the observed SV brings different knowl- 
edge of the core state x compared to the observed field. To see this, let us consider 
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the magnetic induction of the poloidal magnetic field beneath the impenetrable and 
“free-slip” CMB (r = r~) 


b 


m 

l 


-WTT) [Vh -^ B ^ + ri 


dr 2 


1(1 + 1 ) 


(12) 


and in the D” -layer 


br = v d 


dr 2 


i(i + iy 


(13) 


In (12), the subscript “ft” implies the horizontal components of the velocity field v, 
and r] is the magnetic diffusivity of the outer core fluid; rjd in (13) is the magnetic 
diffusivity of the D” -layer (77 < % in general). These two equations show clearly 
that the observed ft" !(0,) will impose the constraint on v, and on the non-potential 
part of the poloidal field. 

The latter, i.e. (13), implies that, at the top of the D”-layer (r = r^), one could 
use a purely potential field ft™^ to match the observed field b]" <0> . However, 
can not recover the observed since 

d 2 bf p) 1(1 + 1 ) m(p) _ 

dr 2 r 2 1 

Therefore, SV assimilation is not redundant to the field assimilation. 

Indeed, our earlier assimilation results in Figure 3 demonstrate clearly that as- 
similation of b r / n< '°' > could not reduce the differences between the forecast SV 
and the observed SV b™^°\ called (O-T) of the SV, although that of the field is 
reduced very rapidly in the first few analysis cycles, a strong indication for the need 
of SV assimilation. 


New Assimilation Approach 

We have been using the sequential assimilation approach in MoSST_DAS (e.g. 
Kuang et al, 2008; Sun and Kuang, 2015). It can be summarized as follows: at the 
analysis time t a when the observation y is made, a new initial condition x a (called 
the “analysis”) is made from the forecast x-f and the observation y, future forecast 
for t > t a can then be made with the following initial value system: 

a-x-f 

— =M(x / ), x / (t a )=x a . (14) 

If there is a linear observation operator H that projects x to the observation space 
(where y is defined), then the analysis x° is of the form 

x“ = x^ + K (y — Hx^) (15) 

K = P^H t (HP / H t + R) 1 (16) 

where K is called the gain matrix, P-^ and R are the error covariances of the forecast 
x-f and of the observation y, respectively. (15) is obtained to minimize the error 
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|H • (x 4 — x a )|“ between the analysis x° and the truth x*. In our previous studies, 
pf is calculated from an ensemble of x-f (Sun et al, 2007; Sun and Kuang, 2015), or 
with some empirical formulations (Kuang et al, 2009; Tangborn and Kuang, 2015). 
The process is repeated again at the next analysis time t a + At (At is called the 
“analysis cycle”). 

If only the observed field is assimilated, then 


y = 


{b7 { °\r d ) 


0 <m<l < L 0 \ = y b . 


By (2) and (5), H is linear and very simple 


(17) 


H = (0,0,H & ,0,0) T , 


(18) 


where corresponds to the subset and has only non-zero entries for b™{rd) 
with 0 < m < l < L a . If the observed SV is also assimilated, then 


y = {yb,y b ) T (19) 

y h = {i>™ (o) M I 0 < m < Z < L 0 } , (20) 

However, by (12) and (13), transformation between and x^ is a differential- 
functional projection and is denoted as H (x^ ) . One could of course construct an 
independent projection system which evaluates H (x-1) directly (e.g. Kalnay, 2003). 
Alternatively, a linearization approximation 'H(x-' ) w H x^ could be made so that 
(15) can still be used. 

There are different means to linearize 'H(x'f). In our current study, we create an 
effective observed field 6™^ defined in the D”-layer that matches both and 

In this approach, comprises of a potential field that accounts for 

and a non-potential field that accounts for 6™^: 


b? M (r) = (^)‘ bf°\r d ) + A- (r - v d f bf°\r d ) 


for r c < r < . (21) 


Obviously, at the top of the D” -layer r = r^, 


,m(o) _ , m(o) 
u l — u l > 


dbf 0) 

dr 


db7 (o) 

dr 


l 


Td 


,m(o) 
U l > 


bt 


= b 


m(o) 

l 


The relative errors of (21) are, via the Taylor expansion, of order [(r<j — r c )/r c ] 3 . 
For a 20km layer thickness, it is smaller than 10 -6 . (21) allows us to extend the 
surface observations to the CMB. The observation vector y is now of the form 


y = 


{b7 {o) {n) | 0 < m < l < L 0 ; r c < n < r d | = y^. 


(22) 


The observation projection is again linear: 

T-L (x^) = H • x-^ 


(23) 
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with H defined in (18). However. Hj, now includes non-zero entries on all grid points 
in the D” -layer r c < r* < r^. 

We can use this approach to further construct an effective observed velocity field 
v° beneath the CMB (r = r~). Since bf 1 is continuous across r = r c the CMB, by 
(12), (13) and (21), we have 


1(1 + 1 ) L 


' h ' 




+ V 


dr 2 


1(1 + 1 ) 


6™ (o) = b, 




(■ r c ) 


= 6™ (o) 


M 




(24) 


Obviously, (24) is an under-determined system, since both b' l nl ' 0 ' > and v° are unknown 
at r~ . But one can find the “best-fit” v° and 6” 1 * 0 - 1 via minimizing the following 
difference 


■»(») 


^jn(o) 

b i (r c ) + 


1(1 + 1 ) 


[V/i • (v h S r )]P - r, 


dr 2 


i(i + iy 


(25) 


If the effective observed velocity field v°(r c ) is included, then the observation vector 

y is 


y = (yv, ys, Yb) T ( 26 ) 

where y^ includes, as shown in (3) and (4), the spectral coefficients c 5™ ° of v° 
at r~. Again, the linearized observation projection (23) is achieved. However, H 
includes additional subsets: 


H=(H„,H*,,H 6 ,0,0) t , (27) 

where H u and H w include only non-zero entries for and at r = r~ , 

respectively. 

Effective Observation Error Covariance 

Since the gain matrix K in (16) depends on the observation error covariance R, we 
need to determine the effective error covariance R for 6™^ which can be calculated 
from those of the Gauss coefficients g™ and /i™. In this section, we only describe a 
formal procedure without going into the details. 

In geomagnetic field modeling (Jackson et al , 2000; Sabaka et al , 2004; Korte and 
Constable, 2005; Olsen et al, 2006; Gillet et al, 2013), the Gauss coefficients, e.g. 
g™, can be described in general as 

9 ™ = S T (t) ■ a lm , (28) 

where S is the vector describing deterministic, model specific base functions in the 
time domain, e.g. B-splinc functions, and a is the coefficient vector which describes 
the observation error statistics. 
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For illustrative purpose, we use the simplest error statistics for our derivation. 
Assume that geomagnetic observations (and thus ot) are unbiased, and with known 
error covariances: 


CX-lm ^ lm 


{e a ) = 0, (e a e^ 


) = C a , 


(29) 


where a\ m = ( ) is the truth (expectation) and C a is the observation error 
covariance matrix of cti m . Thus, by (28), 


g T =g^+e g , 9 r (t) = ST - £g = S T ■ e, 
( e g) = S T ' C a • S = R l ™. 


(30) 


Similar formulation applies to h™ as well. By (10) and (30), we have 

&r (o) 0) = b 7 {t \r) + tb(r), 


b f t \ r) = r l(r i)' Gm 

Cb(r) = y (y) Gm (e g - ie h ) 


This leads to 


( £ fc e b) — 




21 


Gi 


(r 1 ™) 2 + (R i ry 


(31) 

(32) 

(33) 


(34) 


One can use this equation to evaluate the covariance at any location in the mantle, 
including r = the top of the D” -layer. If S in (30) is replaced by S, then we can 
obtain the covariance R l ™ of the SV, 


-qItyi oT P Q 

rig o • C Q • o, 


and therefore the variance of b] 




6P < ”>(r)=6r'(r) + «l(r), 


im(t) , 


( e b e D M = 




21 


Gi 


(■ Rl n +( Rl n 


(35) 

(36) 


The full error covariance of b™^°\r) can then be determined from (21), (34) and 
(36). 


Results 

In this study, we focus only on (22), i.e. assimilation of the effective observed field 
bj 1 ^ which matches both the observed field and the observed SV by^ at the 
top of the D”-layer, mainly for two goals: to explore improvements of the assimila- 
tion system with the observed SV, such as the model spin-up process and rms of the 
observed minus forecast (O-R) of the magnetic field; and to understand responses 
of the core state x to the observed SV, in particular changes of the velocity field v 
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beneath the CMB. Both are critical for determination of the effective velocity field 
v in (25), and thus for implementation of the more comprehensive observation (26). 

We consider only the observations for the time period 1900 — 2000 simply because 
modern observatory and satellite data provide very high quality (g™, h J n ) and (g™, 
hy 1 ). These coefficients are from GUFM1 (Jackson the al, 2000) for 1900-1962 and 
CM4 (Sabaka et al , 2004) for 1962-2000. We also set L a = 8, lower than the highest 
degrees of the two models. For our research purposes, we carry out three distinct 
experiments: 

Case I: Free-running model (no assimilation) 

Case II: Assimilation of with (17) (37) 

Case III: Assimilation of 6™*°^ with (22) 

Except the differences in the data y in analysis, everything else is identical in 
the experiments, including the original initial state at 1900. The analysis cycle is 
At = 5 years. By this design, we can identify exactly the causes of changes in the 
dynamo state x: the differences between the solutions of Case I and Case II are due 
to assimilation of the observed field and the differences between the solutions 

of Case II and Case III are due to the assimilation of the observed SV These 

allow us to understand clearly the responses of the core state to surface observations, 
and their dynamical consequences. 

We use a modeled observation error covariance, since the actual error covariances 
of the field models are not yet available. The model error covariance R is assumed 
diagonal, with the diagonal elements defined as 

Rim = l e fl(0^r| 2 ) £ r(1) = e o{t) + [ £ l(t) — e o(t)} J Ti (38) 

L / 0 -L 

where eo and ei decreases linearly in time: e 0 decreases from 0.01 in 1900 to 0.001 
in 2000, and e\ decreases from 0.3 in 1900 to 0.1 in 2000. These imply that the 
relative errors in (38) decreases in time, but increases with the degree l. 

We would like to point out here that Gillet et al (2013) provided a global field 
model which includes a full error covariance of the Gauss coefficients. This model 
and any future model with specified error statistic knowledge are more appropriate 
for GDAS. However, we conjecture that (38) is sufficient for our current objectives. 

Responses of the Magnetic Field to SV Assimilation 

The quantities used to understand the responses of the magnetic field are the ( O-J 7 ) 
of the radial magnetic field B r and its SV B r . Instead of using traditional (O-T), 
we prefer the following modified definition 
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at r = rd- Replacing by 1 by 6™ in (39), we have (O-T)g of the SV. This modified 
(O-J-) can tell us more accurately how close is the forecast to observation, because it 
eliminates the effect of changes in magnitude of the individual spectral coefficients. 
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Figure 2 are the (0-F) B and (O-J 7 ) B of Case II (dashed lines) and Case III (solid 
lines). From this figure, we can observe clearly that their magnitudes in Case III 
are approximately 30% smaller than those in Case II over the entire assimilation 
period, demonstrating a substantial improvement in forecast accuracies with the 
SV assimilation (21) and (22). 

The SV assimilation also helps accelerate the dynamo model spin-up process. For 
example, we can observe from Figure 2 that the time variations of {0-JF) B are 
nearly identical in both cases: they decay nearly monotonically over much of the 
assimilation period before leveling off in the last 20 years (from 1980 to 2000). But 
(0-T) b , as shown in Figure 3, are very different in the two cases: in Case II, it 
increases first from 1900 to 1940; and only starts to decay continuously in the last 
20 years. In Case III, however, ( O-J 7 ) B decays almost monotonically in time, except 
two small surges around 1940 and 1980. This implies that the dynamo core state 
x-1 responds stronger to the SV assimilation. In other words, the SV assimilation 
helps to accelerate the model spin-up process. 

To better understand how do the forecasts b and 6™^ respond to the obser- 
vations yi, in (17) and in (22), we examine first the ( O-F ) for individual degrees. 
In Figure 4 are ( 0-T) B for the degrees l < 6. Improvements are clearly shown in 
all 6 degrees, as all values are smaller in Case III than in Case II. But we can also 
observe significant differences in individual degrees. For example, ( 0-T) B for the 
odd degrees ( l = 1, 3, 5) increase in magnitude from around 1980. But those for the 
even degrees ( l = 2,4,6) do not show either visible increases or increases far less 
significant than those for the odd degrees. 

As shown in Figure 5, the difference between the odd and even degrees of ( O-J 7 ) B 
is even more significant. There is still a strong surge in magnitude for l = 3 around 
1980 in the both cases. But the reduction for l = 5 is minimal. In particular it 
does not decay monotonically in time in either case. These differences may indicate 
potential inconsistencies between the core dynamics of the model and the time 
variation of the Gauss coefficients. We will discuss this again later in this paper. 

Responses of the Velocity Field to SV Assimilation 

Why does the dynamo model respond faster and stronger in Case III than in Case 
II? We can find at least partial answers from the difference between the free- running 
model solutions x M (Case I), and the forecasts x^ in Cases II and III, in particular 
the differences in the velocity field v beneath the CMB, because they are the direct 
consequences of the magnetic induction (12). The knowledge is also very important 
for obtaining the “effective” observed velocity field (25) for future studies. 

Since in our geodynamo model, the CMB is impenetrable and is free-slip, the 
radial velocity v r = 0 and, by (3) and (4), the horizontal velocity depends on 
dv™ /dr and at r = r c . Therefore, it is very convenient to examine the following 

two variables beneath the CMB: 



(40) 



(41) 
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where v' r is poloidal and describes the up-and-down welling, and uj r is toroidal and 
describes the differential rotation. The rms differences (Ai-A) of these two variables 
between the free-running model solutions (Case I) and the forecasts (Cases II and 
III) can be used to quantify the responses of the core flow to the assimilation of 
surface observations: 


(M-?) vp = Wv^-v^h 
{M-T) vt = \\u,?-u>lh 
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In the above equations, || • || 2 is the L2 — norm (or rms ) over the CMB. 

In Figure 6 are the non-dimensional (with the scaling factor 5 x ICC 6 year -1 for 
dimensional values) \\v’ r \\2 (red) and ||w r ||2 (blue) of the free-running model (Case I). 
As shown in the figure, v' r increases slightly in magnitude in the assimilation period, 
and io r remains flat. But, the rms differences (Ai-A) vp (shown in Figure 7) and 
( AA-A) vt (shown in Figure 8) increase in time, i.e. a growing divergence between 
the forecast state x-' and the free-running model state x M . 

From Figures 7 and 8, we can also observe that (Ai-A) of Case III (the solid lines) 
are slightly larger than those of Case II (dashed lines) , implying that x-^ moves away 
from x M faster with the SV assimilation (22), another demonstration of improved 
model spin-up with the SV assimilation. However, the differences are much less 
significant than those of the magnetic field. This suggests the need for the effective 
observed velocity field v° to increase further (Ai-A) of the velocity field, and thus 
to expedite the model spin-up process. 

To aid the future study of determining the effective core flow from the observed SV 
via (25), we need to understand better the details of ( Ai-A ), e.g. their distributions 
in the spectral space defined by the spherical harmonic degrees l and orders m. 
We shall pay special attention to their distributions in l, i.e. the summation of the 
terms in (42-43) with 0 < m < l for a given degree l; and their distributions in to, 
i.e. the summation of the terms in (42-43) with m < l < Lm for a given order to. 
Since, as shown in Figures 7 and 8, the differences between the two cases are very 
small, we can focus only on Case III without loss of generality. 

In Figure 9 is the distribution of ( A4-A) vp in the degree l, and in Figure 10 is 
its distribution in the order to. From the figures we can find that (A4-A) vp varies 
substantially in the spectral spaces. As shown in Figure 9, the differences for the 
degrees 15 < l < 35 increase the fastest in time, and their magnitudes are the largest 
at the end of the assimilation period, with the peak at l = 20. The differences are 
much smaller and grows slower in time for the degrees l < 5 and l > 40. But, as 
shown in Figure 10, the distribution in to is more broad band: the differences for 
5 < m < 35 increase rapidly in time and reach comparable values in magnitude at 
the end of the assimilation period. However, ( Ai-A ) for to < 4 are very different: 
they remain small and nearly unchanged throughout the entire assimilation. These 
suggest that the responses of the poloidal velocity is dominantly non-axisymmetric 
(to > 0). 
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The distribution of (A / t-J 7 ) VT of the toroidal velocity, as shown in Figures 11 
and 12, displays both similar and distinct characteristics. Its distribution in l is 
very similar to that of (A \-T) vpl except that it peaks at a higher degree l = 30. 
But its distribution in m (Figure 12) is very different: the differences for m < 20 
remain comparable in both the magnitude and the time increasing rate. But they 
decay rapidly for larger m. It should be pointed out in particular that, opposite to 
(M-J 7 ) Vp (in Figure 10), ( M.-F) of the axisymmetric toroidal velocity (m = 0) 

remains very large, implying that the axisymmetric toroidal flow is very sensitive 
to the surface observations. 

Conclusions 

In this study we have examined the consequences of assimilating the observed SV 
on geomagnetic forecasts and on the responses of the dynamo core state. We argued 
that, because geomagnetic data sampling frequencies are several orders of magni- 
tude higher than those of the SV, the geomagnetic field and its SV are concurrently 
measured. We further demonstrated that the observed SV provides unique knowl- 
edge of the magnetic field and the velocity field in the core. Thus assimilations of 
the observed field and of the observed SV are necessary and are not redundant. 

In this study, we incorporate the observed SV into the observation vector y via 
introducing the effective poloidal field (21) in the D”-layer, which is then used 
in the sequential assimilation algorithm (15). We designed three experiments (37) 
to identify the impact of SV assimilation: a free-running model dynamo simulation 
(Case I); an experiment with the assimilation of the observed field (Case II), and an 
experiment with the assimilation of both the observed field and its SV. The relative 
(O-J 7 ) of the field and SV, defined in (39), at the top of the D” -layer are used to 
measure the forecast accuracies; the (M-J 7 ) of the poloidal velocity field (42) and 
of the toroidal velocity field (43) beneath the CMB are used to characterize the 
responses of the core state to the SV assimilation. 

The results of our experiments demonstrate clearly that the SV assimilation with 
(21) improves significantly the geomagnetic forecast accuracies since, as shown in 
Figures 2 and 3, both (0-J 7 ) B and ( O-J 7 ) B in Case III are more than 20% smaller 
than those in Case II. In particular, the improvements occur to all degrees, as 
shown in Figures 4 and 5. The nearly monotonic decay in time of (O-J 7 )^ in Case 
III (Figure 3) shows clearly that the SV assimilation accelerates the spin-up of the 
dynamo model. 

The improvement by the SV assimilation can be also seen from the differences 
(A l-J 7 ) Vp and (A4-A) between the free-running model state and those of the 
assimilations. As shown in Figures 7 and 8, these differences grow rapidly in time, 
showing an accelerated departure of the core state with assimilation from the free- 
running model state. The differences in Case III are slightly larger than those in 
Case II, further demonstrating the improvement brought by the SV assimilation, 
though such increment is less significant that those in the ( O-J 7 ) of the magnetic 
field (Figures 2 and 3). 

Our results have further implications. First, even with the help of (21), the dynamo 
model is still not fully spun up. For example, though (O-J 7 ) B decreases monodically 
in time, the SV forecast is still very far away from the observations, as (0-F) B « 
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0{ 1) for all degrees (see Figure 5). This can be shown further by the continuously 
growing differences ( M-F) vp and (Af-J r )„ T between the forecast velocity field and 
that of the free- running model (see Figures 7 and 8). In addition, the differences 
at the end of the assimilation period are still very small, approximately 10% in 
magnitude of the velocity field of the free- running model (Figure 6). 

These suggest that much larger velocity differences {M-T) vp and ( M-J-) Vt are 
needed over a shorter assimilation period for expediting the model spin-up. Assim- 
ilation of the effective observed velocity (25) could be an answer. However, as dis- 
cussed earlier, (25) is an underdetermined system, since both and (more 
specifically, d 2 b™^°' > /dr 2 ) are unknown beneath the CMB. Thus, the responses of 
v h to the SV assimilation, e.g. (Af-J r ) ?jp (in Figures 9 and 10) and (M-T) vt (in 
Figures 11 and 12) are needed to determine vj°K For example, as shown in the 
two figures, the non-axisymmetric (m > 0) poloidal velocity vj 71 ' around the de- 
gree l = 20, and the toroidal velocity ojJ 71 around the degree l = 30 and order 
m < 20 should be given more attention, as they are most sensitive to the surface 
observations. 

An alternative answer could be the core states inverted from the surface obser- 
vations and dynamo solutions, such as those of Aubert (2013, 2014). These can be 
used as the analysis of the assimilation system. But cautions should be taken with 
this approach. For example, the inverted velocity field beneath the CMB is actually 
derived with the observed field and SV and the magnetic diffusion of the dynamo 
state (Aubert 2014). This could potentially lead to dynamical inconsistencies, as 
well as uncertainties in error statistics. 

Our results also show several new features that may have implications to field 
modeling and to core flow inversion. One new knowledge is from the time variation 
of (0-JV As shown in Figure 5, (O-JF)g of the odd degrees (l = 1,3,5) are 
significantly different from those of the even degrees (l = 2,4,6): the values of 
the even degrees decay nearly monotonically in time; but those of the odd orders 
show either spikes (for l = 1,3) during the assimilation, or even increase over time 
(for 1 = 5). These even-odd degree disparities suggest inconsistencies between the 
model and the observations. These inconsistencies could be entirely due to numerical 
dynamo model which may include a magnetic induction different from those in the 
Earth’s outer core, or may include some mechanisms resulting in different symmetry 
properties of the core state. But the inconsistencies could also come from possible 
biases in the field models that are not included in the observation error covariances. 
For example, ionospheric ring current generated field (an external field component) 
contributes dominantly to the Gauss coefficients of degrees l = 1,3,5, and varies 
on time scales comparable to those of SV, e.g. the solar activity cycles (Sabaka et 
al , 2015). Model biases exist if this part of the signals is not well separated from 
those of the core field. This is potentially an area for application of geomagnetic 
data assimilation. 

The core fluid flow responses, i.e. the differences ( M-T) vp and ( M-F) vt between 
the forecast velocity field v{ and the v M of the free- running model (see Figures 9- 
12), from our experiments could also help inversion of core flow from the observed 
SV. For example, the different characteristics in (M-T) and (M-T) vt suggest 
that the poloidal velocity field and the toroidal velocity field could be treated sepa- 
rately in the core flow inversion. It should be pointed out that purely toroidal core 
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flow approximations were used in previous studies (e.g. Bloxham et a/, 2002; Olsen 
and Mandea, 2008). The strong responses of the high degree core flow (l « 20 for 
the poloidal flow and l ~ 30 for the toroidal flow) to the observed SV (for l < 8) 
indicate that higher degree velocity field should be included in the core flow inver- 
sion. For example, one would normally expect that, due to nonlinear effects, i.e. 
the quadratic terms in Navier-Stokes equation and the induction equation, the core 
flow up to the degrees twice as much as that of the SV should be sufficient for the 
core flow inversion (as in Aubert, 2013). But our results show that time evolution 
of the core flow leads to the strongest responses for the degrees more than triple of 
the maximum degree of the SV. Therefore, inversion of time-dependent core flow 
from the observed SV (up to degree 13) should include high degree (l > 40) spectral 
coefficients. 

Again we should point out that our results could be improved with more sophisti- 
cated assimilation algorithms and field models with more accurate error statistics. 
For example, we anticipate more accurate estimation of (O-T) for both the field 
and the SV, and better assessment of the core state responses if a full ensemble 
approach is used for the covariance P^, and a more appropriate observation error 
covariance, e.g. those determined by Gillet et al (2013), than (38) used in this study. 
Regardless, our assimilation experiments have shown clearly the importance of SV 
assimilation, and the improvements that the SV assimilation brought to forecast 
accuracies and to model spin-up processes. 
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Figures legends 


Figure 1 The time scales tj (1) derived from CM4 for the period 1960-2000: the solid line is for 
the dipole ( l = 1), with 77 « 1500 years; the dashed line is for the non-dipolar components 
(/ > 2), with 77 « 70 years. 


Figure 2 The rms (0-J r ) B of the magnetic field in Case II (dashed line) and Case III (solid line). 
In both cases, (0-.F) s <C 1 and decays monotonically after the first 3 analysis cycles, and then 
levels off in the last 20 years. This shows the continuing improvement in the forecast accuracies. 

In addition, the ( O-J r ) results in Case III (with the assimilation of and 6™^) are in general 

more than 20% smaller than in Case II (with only the assimilation of b™^), showing a clear 
improvement in forecast accuracies. 


Figure 3 Similar to Figure 2, but for ( O-J -) B of the S V. In Case II (dashed line), 

( 0-T) B = 0(1) for much of the assimilation period before decays gradually in the last 20 years, 

implying that there is no similarity between the forecasted SV 6™^ and the observed SV b 
But its magnitude is much smaller in Case III (solid line), and it decays monotonically in time, 
indicating that the SV assimilation accelerates the spin-up process. 
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Figure 4 The of the first 6 spherical harmonic degrees in Case II (dashed lines) and Case 

III (solid lines). 


Figure 5 Similar to Figure 4, but for ( O-J ~) B - 


Figure 6 The non-dimensional ||v £.||2 ( rec Q and ll^vIb/lO (blue) beneath the CMB r = from 
the free-running model solutions (Case I). The dimensional values can be obtained with the 
scaling factor 5 x 10“ 6 year - 1 . 


Figure 7 The (Ai-J 7 ) of the poloidal velocity field v' r as defined in (42). The dashed lines are the 
results without SV assimilation (Case II) and the solid lines are those with the SV assimilation 
(Case III). 


Figure 8 Similar to Figure 7, but for The (A4-J 7 ) of the toroidal velocity u r as defined in (43). 


Figure 9 The distribution of (A4-.F) in spherical harmonic degrees l with the SV assimilation 
(Case III). 


Figure 10 Similar to Figure , but for the distribution of (A4-J 7 ) Vp in spherical harmonic orders m. 


Figure 11 Similar to Figure 10, but for (A / [-J 7 ) Vrp . 


Figure 12 Similar to Figure , but for (A4-J 7 ) Vt . 


